High Reynolds number flow has been simulated by generating discretized solutions of the Navier-Stokes differential equations by performing high-precision floating point arithmetic operations at each of many discrete spatial locations on variables representing the macroscopic physical quantities (e.g., density, temperature, flow velocity). Another approach replaces the differential equations with what is generally known as lattice gas (or cellular) automata, in which the macroscopic-level simulation provided by solving the Navier-Stokes equations is replaced by a microscopic-level model that performs operations on particles moving between sites on a lattice.
The traditional lattice gas simulation assumes a limited number of particles at each lattice site, with the particles being represented by a short vector of bits. Each bit represents a particle moving in a particular direction. For example, one bit in the vector might represent the presence (when set to 1) or absence (when set to 0) of a particle moving along a particular direction. Such a vector might have six bits, with, for example, the values 110000 indicating two particles moving in opposite directions along the X axis, and no particles moving along the Y and Z axes. A set of collision rules governs the behavior of collisions between particles at each site (e.g., a 110000 vector might become a 001100 vector, indicating that a collision between the two particles moving along the X axis produced two particles moving away along the Y axis). The rules are implemented by supplying the state vector to a lookup table, which performs a permutation on the bits (e.g., transforming the 110000 to 001100). Particles are then moved to adjoining sites (e.g., the two particles moving along the Y axis would be moved to neighboring sites to the left and right along the Y axis).
In an improved lattice gas technique, the state vector at each lattice site includes many more bits (e.g., 54 bits for subsonic flow) to provide variation in particle energy and movement direction, and collision rules involving subsets of the full state vector are employed. In a further improved system, more than a single particle was permitted to exist in each momentum state at each lattice site, or voxel (these two terms are used interchangeably throughout this document). For example, in an eight-bit implementation, 0-255 particles could be moving in a particular direction at a particular voxel. The state vector, instead of being a set of bits, was a set of integers (e.g., a set of eight-bit bytes providing integers in the range of 0 to 255), each of which represented the number of particles in a given state.
More recently, Lattice Boltzmann Methods (LBM) use a mesoscopic representation of a fluid to simulate 3D unsteady compressible turbulent flow processes in complex geometries at a deeper level than possible with conventional computational fluid dynamics (“CFD”) approaches. A brief overview of a LBM method is provided below.
Boltzmann-Level Mesoscopic Representation
It is well known in statistical physics that fluid systems can be represented by kinetic equations on the so-called “mesoscopic” level. On this level, the detailed motion of individual particles need not be determined. Instead, properties of a fluid are represented by the particle distribution functions defined using a single particle phase space, f=f(x,v,t), where x is the spatial coordinate while v is the particle velocity coordinate. The typical hydrodynamic quantities, such as mass density, fluid velocity and temperature, are simple moments of the particle distribution function. The dynamics of the particle distribution functions obeys a Boltzmann equation:∂tf+v∇xf+F(x,t)∇yf=C{f},  Eq. (1)where F(x,t) represents an external or self-consistently generated body-force at (x,t). The collision term C represents interactions of particles of various velocities and locations. It is important to stress that, without specifying a particular form for the collision term C, the above Boltzmann equation is applicable to all fluid systems, and not just to the well known situation of rarefied gases (as originally constructed by Boltzmann).
Generally speaking, C includes a complicated multi-dimensional integral of two-point correlation functions. For the purpose of forming a closed system with distribution functions f alone as well as for efficient computational purposes, one of the most convenient and physically consistent forms is the well-known BGK operator. The BGK operator is constructed according to the physical argument that, no matter what the details of the collisions, the distribution function approaches a well-defined local equilibrium given by {feq(x,v,t)} via collisions:
                              C          =                                    -                              1                τ                                      ⁢                          (                              f                -                                  f                  eq                                            )                                      ,                            Eq        .                                  ⁢                  (          2          )                    where the parameter τ represents a characteristic relaxation time to equilibrium via collisions. Dealing with particles (e.g., atoms or molecules) the relaxation time is typically taken as a constant. In a “hybrid” (hydro-kinetic) representation, this relaxation time is a function of hydrodynamic variables like rate of strain, turbulent kinetic energy and others. Thus, a turbulent flow may be represented as a gas of turbulence particles (“eddies”) with the locally determined characteristic properties.
Numerical solution of the Boltzmann-BGK equation has several computational advantages over the solution of the Navier-Stokes equations. First, it may be immediately recognized that there are no complicated nonlinear terms or higher order spatial derivatives in the equation, and thus there is little issue concerning advection instability. At this level of description, the equation is local since there is no need to deal with pressure, which offers considerable advantages for algorithm parallelization. Another desirable feature of the linear advection operator, together with the fact that there is no diffusive operator with second order spatial derivatives, is its ease in realizing physical boundary conditions such as no-slip surface or slip-surface in a way that mimics how particles truly interact with solid surfaces in reality, rather than mathematical conditions for fluid partial differential equations (“PDEs”). One of the direct benefits is that there is no problem handling the movement of the interface on a solid surface, which helps to enable lattice-Boltzmann based simulation software to successfully simulate complex turbulent aerodynamics. In addition, certain physical properties from the boundary, such as finite roughness surfaces, can also be incorporated in the force. Furthermore, the BGK collision operator is purely local, while the calculation of the self-consistent body-force can be accomplished via near-neighbor information only. Consequently, computation of the Boltzmann-BGK equation can be effectively adapted for parallel processing.
Lattice Boltzmann Formulation
Solving the continuum Boltzmann equation represents a significant challenge in that it entails numerical evaluation of an integral-differential equation in position and velocity phase space. A great simplification took place when it was observed that not only the positions but the velocity phase space could be discretized, which resulted in an efficient numerical algorithm for solution of the Boltzmann equation. The hydrodynamic quantities can be written in terms of simple sums that at most depend on nearest neighbor information. Even though historically the formulation of the lattice Boltzmann equation was based on lattice gas models prescribing an evolution of particles on a discrete set of velocities v(ε{ci, i=1, . . . , b}), this equation can be systematically derived from the first principles as a discretization of the continuum Boltzmann equation. As a result, LBE does not suffer from the well-known problems associated with the lattice gas approach. Therefore, instead of dealing with the continuum distribution function in phase space, f(x,v,t), it is only necessary to track a finite set of discrete distributions, fi(x,t) with the subscript labeling the discrete velocity indices. The key advantage of dealing with this kinetic equation instead of a macroscopic description is that the increased phase space of the system is offset by the locality of the problem. Due to symmetry considerations, the set of velocity values are selected in such a way that they form certain lattice structures when spanned in the configuration space. The dynamics of such discrete systems obeys the LBE having the form fi(x+ci,t+1)−fi(x,t)=Ci(x,t), where the collision operator usually takes the BGK form as described above. By proper choices of the equilibrium distribution forms, it can be theoretically shown that the lattice Boltzmann equation gives rise to correct hydrodynamics and thermo-hydrodynamics. That is, the hydrodynamic moments derived from fi(x,t) obey the Navier-Stokes equations in the macroscopic limit. These moments are defined as:
                                                        ρ              ⁡                              (                                  x                  ,                  t                                )                                      =                                          ∑                i                            ⁢                                                f                  i                                ⁡                                  (                                      x                    ,                    t                                    )                                                              ;                ⁢                                  ⁢                                            ρ              ⁢                                                          ⁢                              u                ⁡                                  (                                      x                    ,                    t                                    )                                                      =                                          ∑                i                            ⁢                                                c                  i                                ⁢                                                      f                    i                                    ⁡                                      (                                          x                      ,                      t                                        )                                                                                ;                ⁢                                  ⁢                                            DT              ⁡                              (                                  x                  ,                  t                                )                                      =                                          ∑                i                            ⁢                                                                    (                                                                  c                        i                                            -                      u                                        )                                    2                                ⁢                                                      f                    i                                    ⁡                                      (                                          x                      ,                      t                                        )                                                                                ,                                    Eq        .                                  ⁢                  (          3          )                                    where ρ, u, and T are, respectively, the fluid density, velocity and temperature, and D is the dimension of the discretized velocity space (not at all equal to the physical space dimension).        